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Electromigration-induced flow of islands and voids on the 
Cu(OOl) surface is studied at the atomic scale. The basic drift 
mechanisms are identified using a complete set of energy bar- 
riers for adatom hopping on the Cu(OOl) surface, combined 
with kinetic Monte Carlo simulations. The energy barriers 
are calculated by the embedded atom method, and parame- 
terized using a simple model. The dependence of the flow on 
the temperature, the size of the clusters, and the strength of 
the applied field is obtained. For both islands and voids it 
is found that edge diffusion is the dominant mass-transport 
mechanism. The rate limiting steps are identified. For both 
islands and voids they involve detachment of atoms from cor- 
ners into the adjacent edge. The energy barriers for these 
moves are found to be in good agreement with the activation 
energy for island/void drift obtained from Arrhenius analysis 
of the simulation results. The relevance of the results to other 
FCC(OOl) metal surfaces and their experimental implications 
are discussed. 



I. INTRODUCTION 

Electromigration (EM) describes the biased diffusion 
processes of bulk and surface atoms under the influence 
of an applied electric field (l]-|j| . The driving force is com- 
posed of two major components: one is the direct elec- 
trostatic interaction between the applied field and the 
conductor ions and the second, denoted as the "wind 
force" is due to momentum transfer from the conduction 
electrons impinging upon the ionic cores. The EM force 
F is characterized by an effective charge Z* of the ions, 
which quantifies its strength according to F — eZ*E, 
where E is the electric field. In metals, the wind force is 
the dominant component, thus Z* is negative The 
problem of EM in metal films has been extensively inves- 
tigated experimentally and theoretically for over three 
decades The interest in this problem was moti- 

vated by the fact that EM has been identified as a major 
failure mode of metal interconnects in microelectronic 
devices, where current densities as large as 10 6 A/cm 2 , 
are typical |||l0|]. Indeed, significant experimental effort 
was directed at material and technological issues, aimed 
to reduce EM-induced damage. These issues include the 
choice of metal to be used (e.g. Al, AlCu alloys, or Cu 



dominated alloys), deposition procedures and wire ge- 
ometries. 

Statistical properties of EM phenomena have been 
studied extensively as a function of applied current den- 
sity and temperature. These properties include the time 
to failure of metal wires under current stress Hjqjlll , drift 
velocities jl2],[l3|, resistance- noise @, as well as resis- 
tance increase rates |p| pT| , ^5|Jiq| . Experimental results 
were analyzed using phenomenological equations, such 
as the Black equation ||, that relates the mean time 
to failure t mean to a single activation energy Q through 
tmean oc j~ 2 exp (Q/fc^T), where j is the current density 
and T is the temperature. A similar equation is also used 
to analyze the resistance-increase rate |||l7|,[l8) . This ap- 
proach provides a useful characterization of the wire per- 
formance and reliability. However, it does not yield much 
insight into the fundamental processes that give rise to 
EM. Particularly, the assumption that there is only one 
activation energy is inadequate. This is manifested in 
the dispersion of activation energies reported in the lit- 
erature for a given material. Taking Cu as an important 
example, we find a broad range of reliably measured ac- 
tivation energies, e.g., 0.47 eV ||, 0.7-0.9 eV Jl|, 0.79 
eV @ and 1.21 eV These results indicate that EM 
is not driven by a single atomic diffusion process, but 
is a complex phenomenon that involves a wide spectrum 
of activation energies, and depends on the microscopic 
details of the samples. 

A large variety of EM-related phenomena have been 
observed, such as thinning, void (and hillock) formation, 
followed by their growth and migration [^],[|. The rate 
at which each of these processes develops depends on 
the dominant atomic diffusion mechanism involved. In 
particular, one can distinguish between atom diffusion 
inside crystallites, along grain boundaries, and diffusion 
on surfaces and interfaces. Grain boundary diffusion 
is believed to be dominant in polycrystalline materials 
However, surface (and interface) diffusion becomes 
more important as the cross section of the wire decreases. 
Hence, studies of surface EM at the atomic and nanome- 
ter scales may provide much insight into the fundamental 
EM-induced processes. Recent studies of surface EM on 
single crystal Si samples have focused on EM-induced 
step dynamics [^9|-|lJ. The mobility of steps and dislo- 
cations in thin polycrystalline metal films has also been 
studied [^2 23 1. In theoretical studies using a continuum 
description, the effects of EM on atomic steps |2J] and 
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voids |2^,|2^] were examined. 

In this paper we present a detailed study of diffusion 
of single islands and voids of monoatomic height on the 
Cu(OOl) surface, under electromigration conditions. In 
our analysis we employ a complete set of hopping en- 
ergy barriers for Cu atoms on the Cu(OOl) surface, ob- 
tained using the embedded atom method (EAM) j27j. 
The biased diffusion of islands and voids is studied us- 
ing energy considerations and kinetic Monte-Carlo (MC) 
simulations. The different mechanisms involved in the 
diffusion of islands and voids are identified and the acti- 
vation energies of the corresponding rate-limiting steps 
are quantified. The drift velocities are obtained as a 
function of the EM bias, island/void size and the surface 
temperature. For both islands and voids the slopes of the 
Arrhenius plots obtained from the simulation results are 
found to be in excellent agreement with the activation 
energies of the rate-limiting steps, and thus confirm their 
identification. The drift velocities are found to be linear 
with the bias magnitude. 

The paper is organized as follows. In Sec. II we present 
the model and assumptions. The relevant atomic pro- 
cesses are presented in Sec. III. The simulations and 
results are shown in Sec. IV, followed by a discussion 
and summary in Sec. V. 

II. MODEL AND ASSUMPTIONS 

Surface diffusion is described by thermally activated 
hopping processes of adatoms. The hopping rate h (in 
units of hops per second) of a given atom to each unoc- 
cupied nearest neighbor (NN) site is given by 

h=v cxp(-E B /k B T) (1) 

where v = 10 12 s _1 is the commonly used attempt rate, 
Eb is the activation energy barrier, ks is the Boltzmann 
constant and T is the temperature. 

The activation energy barrier Eb depends on the local 
environment of the hopping atom, namely the configu- 
ration of occupied and unoccupied adjacent sites. We 
assume that only nearest and next-nearest neighbor sites 
have a non-negligible effect on the activation energy. On 
FCC(OOl) surfaces, under these assumptions, the hop- 
ping energy barrier is determined by the occupation of 
seven adjacent sites, as shown in Fig. [I]. In this model, 
adatoms diffuse along the principal crystal directions. On 
the FCC(OOl) surface, the possible paths are along the 
x =< Oil > and y —< Oil > axes. 

In the simulations described below we use a set of en- 
ergy barriers for Cu on Cu(OOl) [^,^9), which includes 
all possible local environments of the hoping atom. The 
effect of the electric current is incorporated by lowering 
the barrier for diffusion in the direction of the electron 
flow, and by raising it in the opposite direction. Thus, 
if the barrier for a certain process is Eb without bias, 
it will be Eb — A for that process in the direction of 



the bias, and Eb + A in the opposite direction. The 
bias A is proportional to the applied field according to 
A = eZ*E ■ a, where a is the projection of the hopping 
distance along the bias direction. In this paper we focus 
on the case in which the field is parallel to either the x 
or the y axis. Thus, for Cu, a — 2.55A, which is the 
lattice constant of the resulting two dimensional square 
lattice. We also assume that bias A is the same for all 
processes and barriers. Although the distortion in the 
energy landscape may depend on the local configuration, 
no details are known about this dependence. It is also 
assumed that the dependence of the attempt frequency 
v on the field is weak enough, and can be neglected. 

Consider a single object (adatom, dimer, vacancy etc.) 
diffusing on the Cu(OOl) surface. The net drift of this 
object in the x direction after time t, will be: Ax = 
n_|_— n_, where n + and n_ are the numbers of hops in the 
positive and negative directions of the x axis respectively. 
On time scales much longer than the time for a single hop, 
the average number of hops will be proportional to the 
elapsed time and to the hopping rate. Assume now that 
the bias is in the positive x direction. According to the 
model we have: 

< n ± >= tv ■ e "b t , (2) 
thus the average drift velocity in the x direction will be: 

<v> =S^ =2 h-smh(^), (3) 

where h = v ■ exp(— EbI^bT} is the rate of the hopping 
process without any bias [Eq. (El)]. The drift velocity of a 
given object is proportional to its hopping rate, namely, 
it depends exponentially on the activation energy, just 
like the diffusion coefficient. The proportionality factor, 
sinh(A/fcsT), is common to all processes, and depends 
on the temperature and the applied field. Under physi- 
cally relevant condition A <C ksT (see below). Therefore 
sinh(A/fcsT) can be approximated by A/fcgT, and the 
drift velocity is linearly proportional to the bias. 



III. ATOMIC DESCRIPTION OF CLUSTER 
DRIFT 

A. The Energy Barriers 

Within the assumptions described above, adatom dif- 
fusion on the Cu(OOl) surface involves 2 7 = 128 hop- 
ping processes (Fig. pi). In order to obtain the rates of 
these processes one needs to calculate their energy bar- 
riers. Our analysis is based on semi-empirical calcula- 
tions obtained via the embedded atom method (EAM) 
This method provides a good description of self 
diffusion of Cu on Cu(OO l) pql , as well as self diffusion 
on other metal surfaces [pO| , pif . Specifically, we use the 
EAM functions of Cu developed by Adams, Foiles, and 
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Wolfer |32| which are fitted to a similar data base as the 
one employed by Foiles, Baskes, and Daw J3^]. Gener- 
ally, these energy barriers are in agreement with those 
obtained using other semi-empirical methods such as the 
effective medium theory (EMT) |34|,|3j| , and ab-initio cal- 
culations Though each process has a different ac- 
tivation energy, they can be roughly divided into four 
These groups give rise to four typical 



groups |29|]37 

time scales in the problem which span over several orders 
of magnitude. These time scales correspond (in ascend- 
ing order) to attachment processes, edge diffusion, free 
surface hopping and detachment processes. We will ar- 
gue that all four scales are essential in the description of 
cluster drift. Furthermore, even much smaller differences 
between processes which belong to the same group will 
appear to be significant. 



B. Simple Model for Energy Barriers 

The calculated energy barriers can be well approxi- 
mated by a simple model which has four parameters [ p9[ . 
This model provides a systematic description of the en- 
ergetics of the hopping processes, and a classification of 
these processes to four groups. It also enables generaliza- 
tion to different materials. According to the model, the 
energy barrier Eg for a certain process depends on the 
seven sites of Fig. [j] as follows: 

Eb = Eq + E^ N ■ (S3 + Si + S5) 

~^NN ' 0^1 + ^5 + 5*2 + 5*6) 

+E NNN ■ (So + S 2 + Si + S 6 ) (4) 

where Si = 1 if site i is occupied and Si = if it is 
vacant, Eq is the barrier for an isolated adatom hopping 
on the surface, E™ N and E^nn are the effective binding 
energies of the hopping atom to its nearest and next- 
nearest neighbors, respectively. The highest energy along 
the hopping path is obtained in the vicinity of the bridge 
site. The binding energy of the hopping adatom at this 
highest point, to an atom in one of the four adjacent sites 
is given by E l ^ p N . Here we use a simplified version of the 
model assuming that E™ N ~ E N y N . This approximate 
degeneracy is found for most FCC metals [£9[|37|]. The 
expression we use for the barriers is thus 
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B = Eq + Enn ■ (S3 



-Ennn ■ (So 



— S2 — Se) 
5*2 + S4 + Sq). 



(5) 



The three parameters Eq, Enn and -Eawat were esti- 
mated for Cu(OOI) and yielded E = 0.49 cV, E NN = 
0.27 eV and E NNN = 0.027 eV H& While sim P le and 
intuitive, the barriers obtained from the model may de- 
viate in some cases significantly from the EAM barriers. 
Such deviations occur in relatively dense local environ- 
ments, where several sites in both sides of the hopping 
atom are occupied. The model accounts for the complex 



interactions between these atoms only on the average. In 
the following discussion barriers obtained by both EAM 
and the model of Eq. (||) are quoted. The data for the 
simulation results presented below were obtained using 
the EAM barriers. 



C. Physical Conditions 

The following discussion concentrates on clusters of 60- 
1000 atoms/vacancies on the Cu(001) surface. Such clus- 
ters are typically created during deposition and sputter- 
ing experiments |3^ , |39|] . We consider temperatures in 
the range 200-600K, which is the relevant range for most 
experimental studies. The model for diffusion described 
above applies well in this regime, where no other hop- 
ping mechanisms are significant. At higher temperatures 
additional processes may take place [ p6[ . 

The effective charge for EM in bulk copper was found 
experimentally to be Z* w —5 |Q. For surface EM 
there exist only theoretical estimations. These calcula- 
tions found the effective charge to be of the same order 
of magnitude, Z* « —20 [§,0, almost independent of 
surface orientation and diffusion path. Typical current 
densities in EM experiments (and in integrated circuits) 
are j « 10 6 — I0 7 A/cm . The resistivity of Cu under 
typical conditions is p w 2 ■ 10 -6 f2-cm. If we assume that 
the current is distributed uniformly across the conduc- 
tor's cross-section, we find the EM force acting on surface 
atom to be F = eZ*jp « 10 eV/cm. For Cu(001) the 
distance between adjacent sites is w 2.55A. The work 
done by the EM force during one hop in the bias di- 
rection is thus w 10 _7 eV. Intuitively, this would be the 
change in the energy barrier for the process. Since most 
of the relevant barriers are in the range 0.1-0.8 eV, we 
find that A is typically 6-7 orders of magnitude smaller 
than Eb- For such a small bias to produce any effect an 
extremely long time is required. In experiments, it takes 
at least several hours at temperatures higher than 400° C 
for any significant change. Simulating systems at these 
temperatures for more than a second is beyond the power 
of contemporary computers. To overcome this difficulty, 
we use bias values in the range A w 10~ 3 — 10~ 4 eV. 
These values are larger than the realistic ones, but still 
much smaller (2-3 orders of magnitude) than the diffu- 
sion barriers. The linear response approximation is thus 
still valid. 



D. Island Drift 

During the drift, clusters (islands and voids) maintain 
an approximate square shape with small fluctuations and 
rounded corners. There is experimental evidence that 
this is the equilibrium shape of islands ]38|j3l| . These 
experiments show that even if a different shape is cre- 
ated (during coalescence, for example), it rearranges to 
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the square-like pattern within several minutes at room 
temperature. This is a consequence of the square sym- 
metry of the lattice and the fast edge diffusion. This 
fact has important implication on the atomic details of 
island diffusion, with or without bias. When clusters dif- 
fuse, they typically move one lattice site at a time, while 
maintaining their equilibrium shape in the new position. 
Although the bias drives the system out of equilibrium, 
since the applied bias is small, the basic pattern is pre- 
served. This feature is also observed in our computer 
simulations. However, in the presence of bias, atoms drift 
along island edges with a preferable direction, and thus 
sharpen the rounded corners in the bias direction. The 
basic cycle of island drift, in which the center of mass 
moves one lattice site in the bias direction, can be divided 
to three main stages. These are shown schematically in 
Fig. 0. To describe the basic cycle, take the starting 
point to be an island with a straight edge in the direc- 
tion of the drift. [Fig. ||(a)]. First, an atom detaches 
from one of the four corners and arrives at the island 
front [Fig. 0(b)]. Next, several other atoms nucleate to 
this atom from the corners and create a new row of atoms 
in the front [Fig. §(c)]. Last, atoms from the rear side of 
the island fill the place of the atoms that formed the new 
row and straighten the corners to form a configuration 
similar to that of Fig. |(a), shifted by one lattice site in 
the bias direction. The entire process can now repeat as 
the drift goes on. We will now go into the details of each 
stage. 

Consider an island with a straight front edge in the 
drift direction as in Fig. ||(a). One atom now detaches 
from one of the corners and arrives at this front. If it 
detached from a rear corner (with respect to the drift 
direction), it will hop along the sides parallel to the drift 
under the bias influence, and arrive at the front. There 
are four possible moves in which an atom can leave the 
corners, which we index from (a) to (d), as shown in 
Figs. |^(a-d), respectively. All these processes actually 
take place, but we will now argue that one of them is 
much more probable, and thus determines the properties 
of the drift. This process is an escape from a straight 
corner [Fig. ||(a)]. By escape here we mean the detach- 
ment of an atom from a relatively bound configuration 
to the island edge, but not a complete detachment from 
the island. 

The activation energy for this move is given by: 

E eS cape(a) = E + Enn + EjyNN 

= 0.79cV(model),0.78 eV (EAM), (6) 

where the first number corresponds to the model [Eq. 
(||)], and the second number is the EAM barrier. After 
the first of the two hops which compose the move in Fig. 
||(a) has taken place, there are even chances (neglecting 
the bias) for that atom to go back to the corner or to 
actually escape to the straight edge. Another possibility 
is the escape from a straight edge [Fig. ||(b)] . The barrier 
for such event is 



E eS cape(b) = E + £j\TJV + ZE^NN 

= 0.82eV(model),0.90 eV (EAM). (7) 

It may seem at first sight that the extra NNN bond in 
comparison to the corner escape may be compensated 
by the large number of atoms at this configuration along 
the island edges. This is not the case, however, because 
an actual escape can happen only one lattice position off 
the corners, and not anywhere along the edges . Oth- 
erwise, the escaped atom is most likely to return to its 
initial position. Furthermore, as mentioned before, the 
barriers obtained by the model are less accurate at dense 
environments. The corresponding EAM barrier for the 
escape from straight edge is 0.9 eV. Another possibility 
is the escape from a rounded (kinked) corner [Fig. ^(c)]. 
The expression obtained from the model for this barrier is 
the same as in Eq. (Q) . The EAM barrier for this move is 
E escape (c) = 0.83 eV, which is slightly higher than that 
for move (a). Also, at this configuration several other 
processes are much more likely to happen before an es- 
cape move takes place. The last move (d) is an escape of 
an atom from an edge kink followed by hopping around a 
corner [Fig. ||(d)]. The barrier for each of these processes 
separately (Ek and E c respectively) is lower relative to 
the previous three: 

Ek = E + 2Ennn 

= 0.54 eV(model),0.48 eV(EAM), (8) 

and 

E c = E + Ennn 

= 0.52 cV(model),0.54 eV(EAM), (9) 

To understand why this move is relatively unlikely, it 
should be noticed that in most cases the atom will re- 
attach to the kink rather than hop around the corner. 
The barrier for a move from the corner back towards 
the kink site is low: Eb ac k = Eq — Enn + Ennn [=0.24 
eV(model),0.18 eV(EAM)]. The average number of times 
the atom has to detach from the kink until it hops once 
around the corner is given by 

(h c + h back )/h c ~ exp \{E C - E back )/k B T], (10) 

since h c <C h ba ck for the relevant temperature range. The 
prefactor depends on the distance between the kink and 
the corner. The effective barrier for the move shown in 
Fig. |(d) is 

E e scape(d) = Ek + E c — Eback 

= E + Enn + 2Ennn 

= 0.82eV(model),0.84 eV (EAM). (11) 

The next stage in the island drift process is the nu- 
cleation of other atoms to the atom hopping along the 
straight facet. The basic mechanism for this nucleation 
is shown in Fig. 0. The energy barrier for this process is 
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Enuc = Eq + 3EnnN 

= 0.57 eV(model),0.53 eV(EAM), (12) 

This process must occur several times in order to create 
a stable new row of atoms on the existing facet. Still, the 
time scale of this nucleation process is short compared to 
that of the escape process. The ratio between the rates 
of these processes is at least 

exp [(E eS cape - E nuc )/k B T] ~ 

~ exp [{E NN - 2E NNN )/k B T}. (13) 

For Cu, this ratio is 125 at 600 K and more than 15000 
at room temperature. In the clusters we consider, where 
a typical linear island size is of 10-30 atoms, this is a 
relatively fast process. 

Now, atoms from the edges parallel to the bias direc- 
tion will fill-in the kinks created by the atoms that moved 
to the new row. This stage is necessary to enable the nu- 
cleation of the next new row. This is due to the high rate 
of adatom re-attachment to kinks. There are two basic 
mechanisms for this edge drift which are shown in Fig. 
pi. In Fig. ||(a) a single atom hops along the edge in the 
bias direction. Eventually this atom will be absorbed in 
a kink site at the island corner. In Fig. |J(b) the same 
happens with an edge vacancy in the opposite direction. 
The barriers for these processes are 

E a = E — Ejvn + 2Ennn 

= 0.27 cV(modcl), 0.25 cV(EAM), (14) 

and 

Eb = Eo + 2Ennn 

= 0.54 cV(model),0.48 eV(EAM), (15) 

The first one is of course much faster, but requires an 
escape move to begin with. To estimate the time scale of 
these processes, one needs to obtain the average number 
of single hopping events required for one atom/vacancy 
that leaves one corner to be absorbed at the opposite 
one. It is reasonable to treat both atoms and vacancies 
as one dimensional random walkers. The problem then 
reduces to a random walk with two absorbing barriers. 
The expression obtained from the theory of biased ran- 
dom walks (see e.g. Ref. [fy]), depends on the bias and 
the length of the walk (i.e. the linear size of the island). 
In the range of these two parameters we deal with, the 
results are 26-200 events. Again, this is a short time scale 
compared to that of escape events. In summary, the drift 
of islands can be divided into three stages with three dif- 
ferent time scales. The velocity of the drift is determined 
by the slowest of these, namely the escape event. We 
come to the conclusion that the escape process (a) is the 
rate limiting step of island drift. The effective barrier for 
the drift is thus 

E B (island drift) = E escape (a) 
= E + Enn + Ennn 
= 0.79eV(model),0.78 eV (EAM). (16) 



E. Void Drift 

Basically, the three stages that exist in the island drift 
process, are found for vacancy clusters as well. The 
atomic details however, reveal important differences. The 
most essential difference is in the equivalent process to 
the corner escape move. Consider a void with a straight 
corner (Fig. ^). The barrier for this detachment process 
is 

Edetach = E + 3EnnN 

= 0.57 cV(modcl),0.53 eV(EAM). (17) 

In most cases however, the detached atom does not es- 
cape from the corner, but goes back to its initial position. 
There are several possible moves in the nearby environ- 
ment that may prevent the detached atom from moving 
back. The one with the lowest barrier is the replacement 
move shown in Fig. ^[ The average number of times a 
detachment event should occur before such move takes 
place, is very well approximated by the ratio of the rates 

hback I h-replace = exp [(E rep i ace — E oac k) / k B T]. (18) 

The effective barrier obtained for an escape from the cor- 
ner is 

E 'escape Edetach E re place Eback 

= E + E NN + 3E NNN 

= 0.84eV(model),0.73 eV (EAM). (19) 

The nucleation of other vacancies on the void edge is 
created mainly by the detachment of atoms from an edge 
kink and their diffusion along the edge (Fig. |7|). The 
bottleneck of this process is the move across the corner 
of the void, shown in Fig. f?](b). The barrier for leaving 
the corner is E c = E + 3Ennn [=0.57 eV(model),0.53 
eV(EAM)]. 

E c = E + SEnnn 

= 0.57 cV(modcl),0.53 cV(EAM), (20) 

In contrast to islands, void corners are "attractive". 
Since the barrier for leaving the corner is comparable 
to the barrier for leaving edge kink, in many cases there 
will be no accumulation of atoms near the corner, and we 
will not get the rounded corners as in islands. In other 
words, the stages of nucleation and corner straightening 
are combined in voids. In order to compare the time 
scales of the stages, we need to estimate how many cor- 
ner detachment events are required for one atom to reach 
the opposite edge of the void. (The hopping on the edge 
itself is again very fast and its actual time is negligible.) 
Using again the approximation of a one dimensional ran- 
dom walk, we obtain an estimation of « 10 such events 
for voids of linear size 10. Since all atoms on the facet 
must go through the corner, at least 100 events are re- 
quired for the nucleation of a complete row. The time 
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scale of the nucleation stage is thus short compared to 
that of the escape event at low temperatures. Under 
these conditions the effective barrier for void drift is 

.Eg (void drift) = E escape = Eq + Enn + 3Ennn 

= 0.84eV(model),0.73 eV (EAM). (21) 

At higher temperatures, however, depending on the size 
of the cluster and the bias, the two time scales may be- 
come comparable. 



F. The Drift Velocity 

For both islands and voids we identified the rate lim- 
iting processes and their energy barriers. These barriers 
are the effective activation energies of the entire drift pro- 
cess. We thus expect the drift velocities of islands and 
voids to scale like 



Vdrift (island) ~ exp(- 



Eq + Enn + Ennn 
k B T 



) (22) 



and 



Vdrift (void) ~ exp(- 



Eo + Enn + 3Ennn 
knT 



) (23) 



respectively. It is convenient to measure the velocity of 
the cluster center of mass in units of lattice sites per 
second. We expect from Eq. (||) that the drift veloc- 
ity is linearly proportional to the bias magnitude. It is 
also proportional to the probability of an escaped atom 
to form a new row rather than go back. This probabil- 
ity cannot be deduced from simple arguments, since it 
generally involves a random walk with moving bound- 
aries and many particles. In general it may depend on 
the cluster size and the temperature. For larger clusters, 
more atoms need to nucleate to complete a new front 
row. It also takes longer to fill the kink site created by 
each newly nucleated atom, since the sides parallel to 
the bias component are now longer. We thus expect that 
the prefactor decreases as the cluster size increases, but 
cannot provide a specific functional form. 



IV. SIMULATIONS AND RESULTS 

In the simulations reported below we used the con- 
tinuous time kinetic MC technique |42]-[48|]. This tech- 
nique is particularly suitable for the simulation of non- 
equilibrium processes, keeping track of the physical time 
in a realistic manner. During the kinetic MC simula- 
tion, the next move is selected randomly from the list of 
all possible moves at the given time with the appropriate 
weights. The time is advanced after each move according 
to the inverse of the sum of all rates. All the processes 
allowed by the model are incorporated in the simulations 
with the appropriate energy barriers. In the simulation 



results presented below we used barriers for Cu/Cu(001) 
obtained by EAM §§|. 

We have performed systematic simulations on single 
isolated islands and voids. The initial island and void 
configurations were chosen to be of a square shape from 
the reasons already discussed. The morphology and lo- 
cation of the clusters were followed for different temper- 
atures, bias directions, and bias magnitudes. This was 
repeated for different sizes of the islands/ voids. As ex- 
pected, islands drift in the direction of the bias, while 
vacancy clusters drift in the opposite direction. In the 
temperature range considered here (200 — 500K), both 
islands and voids drift as a whole, since the activation 
energy for a complete detachment of atoms or vacancies 
is high relative to the processes discussed above. Fur- 
thermore, even if an atom/vacancy is detached from the 
cluster, the bias is not strong enough for a substantial 
drift, and it will re-attach after few moves. 

In the simulations we follow the displacement of the 
center of mass of the clusters as a function of the physical 
time. It was found that under the EM bias, clusters drift 
on the average at a constant velocity. The dependence of 
this velocity on the cluster size, bias magnitude and the 
temperature is shown below. 



A. The Effect of Bias Magnitude on the Drift 
Velocity. 

We found that the drift velocity of islands depends 
linearly on the bias as can be seen from Fig. ||. This de- 
pendence is expected from Eq. (|j|) , and was confirmed by 
simulations for bias magnitudes in the range A m 10~ 3 - 
10~ 5 cV. In addition we confirmed by simulations that 
the effect of bias direction is given by v x oc A cos(</>) and 
v y cx Asin(^), where v x and v y are the x and y compo- 
nents of the drift velocity, and A is the bias magnitude 
applied in an angle (f> relative to the x axis. 



B. Drift Velocity as a Function of the Temperature 

The dependence of the drift velocity on the tempera- 
ture is shown in Fig. ^. The activation energy for island 
diffusion, is found by the best exponential fit to be 
£ B (island drift)= 0.78 ± 0.02 eV. This coincides with 
the energy barrier for the escape move, which we iden- 
tified as the rate limiting step in the analysis above. 
For voids, we find that for temperatures in the range 
22QK < T < 300K there is a good fit to the predicted 
value: 

i?e(void drift)= 0.73 ± 0.02. At higher temperatures 
there is a small deviation from this value towards a lower 
value. This may indicate that the probability of an atom 
that escaped from a void corner to be re-attached rather 
than start a new vacancy row, slightly depends on the 
temperature. 



G 



C. Dependence of The Drift Velocity on Cluster Size 

Fig. shows the dependence of the drift velocity on 
the linear size of the cluster. As could be expected, there 
is a monotonic decrease in drift velocity. The quantita- 
tive details, however, are not completely concluded. 

V. DISCUSSION AND SUMMARY 

In this paper we have studied the mechanisms of EM- 
driven diffusion of single islands and voids on the Cu(OOl) 
surface. We found that the drift velocity of a cluster 
of a given linear size depends on the bias A and the 
temperature T according to 

v(T,A)=A--^-vexp(-Q/k B T), (24) 
k B l 

where Q is the activation energy of the rate limiting step 
of the drift process and A is a constant. The reciprocal 
dependence on the temperature could not be inspected 
from the simulation results, since it is much weaker than 
the exponential dependence. It is deduced, however, from 
the discussion following Eq. (0), and from dimensional 
analysis. The activation energy Q is given by Eq. 
for islands and by Eq. (|2l]) for voids. 

The value of the constant A is proportional to the prob- 
ability of an escaped atom to form solid nucleation and 
establish a new front row of atoms. As we mentioned be- 
fore, the atomic diffusion processes that determine this 
probability are complicated. The simulations, however, 
yield a value of A ~ 5 — 20 for islands, and A ~ 40 — 80 
for voids, depending on the cluster size, provided that the 
velocity is measured in units of lattice sites per second. 

It was found |37j that most FCC(OOl) metal surfaces 
share the same qualitative features of the different hop- 
ping mechanisms. The hopping energy barriers in these 
systems can be parameterized by models similar to the 
one we used here, with the specific parameters appropri- 
ate for each metal. The analysis of the atomic processes 
involved in the drift of islands and voids, and their typ- 
ical time scales should be thus valid for these FCC(OOl) 
metal surfaces as well. Several simulations performed 
for Ag/Ag(001) diffusion, indicate that this is indeed the 
case. 

Extrapolation of the results to other regimes are rea- 
sonable as long as the relations between the time scales 
remain as in the above discussion. Larger clusters or 
weaker bias may be considered as long as the nucleation 
and corner straightening stages are still faster relative to 
the escape move. For islands at room temperatures, for 
example, this means clusters of up to ~ 10 4 atoms at bias 
of 0.001 eV, or bias down to ~ 10~ 6 eV for clusters of 
~ 100 atoms. The linear dependence of the drift veloc- 
ity on the bias makes the extrapolation down to a more 
realistic bias values straightforward. 



In order to test our predictions experimentally one 
needs to prepare Cu(001) surfaces with islands of the 
desired size distribution. Then one needs to drive cur- 
rent along the surface at high current densities. Using 
scanning tunneling microscope (STM) one will be able 
to measure the drift velocities of different islands and 
obtain their dependence on the island size, bias and tem- 
perature. Related experimental work has been done on 
diffusion of islands on Cu(001) with no bias and the 
dependence of the diffusion coefficient on island size was 
obtained. To our knowledge, no such experiments on 
single crystal metal surfaces have been done under EM 
conditions. It seems that a serious difficulty in perform- 
ing such experiments on single crystal samples may be 
to reach the high current densities required to obtain a 
fast enough drift, due to the considerable width of such 
samples. 
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FIG. 1. Classification of all possible local environments of 
a hopping atom including seven adjacent sites. Each site can 
be either occupied or unoccupied, giving rise to 2 7 = 128 local 
environments. Sites 1, 3 and 5 are nearest neighbors of the 
original site while sites 1, 2, 5 and 6 are adjacent to the bridge 
site that the atom has to pass. 

FIG. 2. The main stages in cluster drift starting from a 
straight facet (a). An atom is detached from the corner and 
starts a new row (b) . Other atoms nucleate to the first atom 
and complete the row (c). Atoms then drift along the edges 
and retain the straight facet . The arrow on the grid indicates 
the initial position of the lattice front. 

FIG. 3. The main mechanisms for detachment that starts 
a new row: (a) from a straight corner, (b) from a kinked 
corner, (c) from a straight edge and (d) from a kinked edge. 
It is found that (a) has the lowest activation energy and is 
thus dominant. 

FIG. 4. The main mechanism for nucleation of atoms on a 
new row. Notice that it can take place only near the corners, 
and that the inverse process is identical. 

FIG. 5. The main mechanisms for atom drift along island 
edges, (a) A single atom hops on a straight edge; and (b) 
an edge vacancy drifts in the opposite direction. The arrow 
indicates the bias direction. 

FIG. 6. The atomic process that starts a new vacancy row 
in a void drift. An atom is detached from the straight corner 
. Another atom may fill its place and create a stable dimer. 

FIG. 7. Nucleation of vacancies at the new vacancy row 
is actually composed of atom detachment from an edge kink 
(a), hopping to the corner of the void (b), and then to the 
opposite corner. 

FIG. 8. Drift velocity as a function of bias magnitude for 
10 x 10 clusters. Stars are voids at temperature of 350K and 
Circles are islands at 400K. Error bars are about the size of 
the symbols. 

FIG. 9. Drift velocity as a function of T _1 on 
semi-logarithmic scale, for 10 x 10 island (stars) and void (cir- 
cle) with bias of 0.003 eV. Solid lines are a fit to an exponential 
form A ■ exp(£^). The temperature range is 300-400K and 
the error bars are smaller than the symbols. The slopes ob- 
tained from the fit are Q = 0.78 ± 2 eV for islands. For voids 
at T < 300 K Q = 0.73 ±2 eV. These values are in agreement 
with the barriers of the rate limiting steps discussed in the 
text. 



FIG. 10. Drift velocity as a function of cluster size for (a) 
islands at bias of 0.002 eV and (b) voids at bias of 0.0005 eV. 
The temperature is 400K. The error bars in (a) correspond 
to one standard deviation of the runs, while in (b) they are 
smaller than the symbols. 
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